A1 = readmatrix('NiI2_bulk_20Vcool_0Vwarm.txt')
A2 = readmatrix('NiI2_bulk_-20Vcool_0Vwarm.txt')

x1 = A1(:,1)
y1 = A1(:,2)
m1(1) = x1(1)
n1(1) = 42

% for temperature below 55K, we binned the data point in the integration.
% coefficient 20 is the sample dimension 2s/(1mm*0.1mm). Size is 1mm here beacuse
% the edge of sample was cut to expose the fresh surface.
for i=2:9
    for j=1:10
        m1(i)=x1(i*10-9)
        n1(i)=(y1(i*10-9+j)+0.005)*20+n1(i-1)
    end
end

for i=10:40
    m1(i)=x1(i+80)
    n1(i)=(y1(i+80)-0.023)*20+n1(i-1)
end

for i= 41:43
    for j=1:10
        m1(i)=x1((i-28)*10-9)
        n1(i)=(y1((i-28)*10-9+j)-0.02)*20+n1(i-1)
    end
end

% coefficient 15 is the sample dimension 1.5s/(1mm*0.1mm)

x2 = A2(:,1)
y2 = A2(:,2)-0.02

m2(1) = x2(1)
n2(1) = -44
for i=2:9
    for j=1:10
        m2(i)=x2(i*10-9)
        n2(i)=(y2(i*10-9+j)-0.01)*15+n2(i-1)
    end
end

for i=10:40
    m2(i)=x2(i+80)
    n2(i)=(y2(i+80)+0.005)*15+n2(i-1)
end

for i= 41:43
    for j=1:10
        m2(i)=x2((i-28)*10-9)
        n2(i)=(y2((i-28)*10-9+j)+0.01)*15+n2(i-1)
    end
end

hold on

yyaxis right
fill(x1,y1,[0.8,0,0],'EdgeColor','none','FaceAlpha','0.7')
fill(x2,y2,[0,0,0.8],'EdgeColor','none','FaceAlpha','0.7')
xlim([50, 65])
xlabel('Temperature [K]')
ylabel('Pyroelectric Current [pA]')
set(gca,'XMinorTick','on','YMinorTick','on')
ylim([-1.3,1.3])
yticks([-1:0.5:1])

yyaxis left
plot(m1,n1,'-','Linewidth',1,'Color',[0.8 0 0])
plot(m2,n2,'-','Linewidth',1,'Color',[0 0 0.8])
ylabel('Polarization [uC/m]')
ylim([-60,60])
yticks([-40:20:40])
xlim([50, 65])
set(gca,'XMinorTick','on','YMinorTick','on')

ax = gca;
ax.YAxis(1).Color = 'k';
ax.YAxis(2).Color = 'k';

set(gca,'box','on','Xcolor',[0 0 0],'Ycolor',[0 0 0])
legend('20 kV/m', '-20 kV/m')
legend boxoff
hold off

set(gcf,'Position',[0, 0, 230, 200])